This was used to determine the specular reflectance and global absorption. 



clc; 
source = "source";
theta0 = 35; # angle of reflectance (negative angle of incidence)
Theta = 1.2526 * 3.14159 / 180; # for Pump

f = 3; #frequency point

####################################################################
# far-field setup                                                  #
####################################################################

## get far-field data
m="reflectance"; # name of monitor measuring near-field reflectance
res = 1000; # x-y resolution of far-fiel
far_E2 = farfield2d(m,f, res); # get far-field with frequency point f
theta = farfieldangle(m,f, res); # get ux
far_total = farfield2dintegrate(far_E2,theta); # get integrated total |E|^2 in far field
plot(theta,far_E2);

# beam parameter
w0 = 15 * 1e-4; # from um to cm
z_R = w0 / Theta; # Rayleigh length in cm
r = 7.5; # in cm, distance of circular detector from sample (if sample surface is at z=0, else recalc)
w_10 = w0 * sqrt(1 + (r / z_R)^2); # in cm is beam radius at 10 cm distance (theoretical)
D = 1.2; # diameter of first mirror in cm

# expected spot position (as incidence)
theta_var=atan(D/(2*r)) * 180 / pi; # far-field angle for variation around theta0 (+-)
#theta_var = 26.74/2* 180 / pi; 
# calculate integrated far-field within solid-angle of detector
# has to be understood that it is "rotation" symmetric from 0 to D
far_R_spec = farfield2dintegrate(far_E2,theta,theta_var,theta0)/far_total;
far_R_diff = 1 - far_R_spec;
?far_R_spec;

R = getresult("reflectance", "T");
R = R.f; 
#R = 1-T.T(f); # Total reflection

real_R = R * far_R_spec; #comparable to experiment


############################
# Determine absolute absorbed Power 
############################
################################################
# calculate power absorbed in  region where  
# the refractive indes is that of a particular materials
# for all frequencies,
# by integrating the absorption over this material.
#################################################

# get monitor data
runanalysis;
m="pabs_adv";
Pabs=getresult(m,"Pabs");
n = getresult(m+"::index","index");
n = n.index_x;
f2=Pabs.f; #get frequencies 
nx=length(Pabs.x);
ny=length(Pabs.y);
nz=length(Pabs.z);
nf=length(Pabs.f);

E = getresult(m+"::field", "E"); 
E2 = E.E2; 
x = E.x; 
y = E.y; 

# define material filters

filter = "Al 1040nm";
filter_real = matrix(nx,ny,nz,nf);
filter_imag = matrix(nx,ny,nz,nf);

n_filter_real = real(getfdtdindex(filter,f2,min(f2),max(f2)));
n_filter_imag = imag(getfdtdindex(filter,f2,min(f2),max(f2)));

#relative difference for almosequal comparator
rel_diff = 1e-15;

for (i=1:nf) {
    filter_real(1:nx,1:ny,1:nz,i) =  almostequal(n_filter_real(2),pinch(real(n),4,2),rel_diff);
    filter_imag(1:nx,1:ny,1:nz,i) =  almostequal(n_filter_imag(2),pinch(imag(n),4,2),rel_diff);
}
filter=filter_real*filter_imag; #matrices zero for vacuum

E2sur = E2*filter;

#integrate
Pabs = integrate2(Pabs.Pabs*filter,1:3,Pabs.x,Pabs.y,Pabs.z);